##########################
#Alternative Explanations#
##########################
#This to obtain:
#table_A4.tex

library("broom")
library("dplyr")
library("grid")
library("tidyr")
library("sandwich")
library("lfe")
library("broom")
library("dplyr")
library("tidyr")
library("grid")
library("stargazer")
library("lmtest")
library("multiwayvcov")
library("ggplot2")
library("glue")
library("spdep")
library("spatialreg")
library("stringr")
library("readxl")  
library("rgdal")


setwd("/Users/bgpopescu/")
#setwd("C:/Users/bogdanp/")


############
#1857 Data#
###########

census_1857 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                       layer="census_1857_points")

census_1857<-subset(census_1857, country=="Croatia" & (treat == 1 | treat == 0))
census_1857$treat<-0
census_1857$treat<-ifelse(census_1857$kreis_comitat_regiment_latin == 'Militar-Grenz-Communitaten' | census_1857$kreis_comitat_regiment_latin == 'Militar-Grenz-Regiments-Bezirke', 1, 0)

census_1857$dist1 <- as.numeric( with (census_1857,ifelse(census_1857$treat==1, 1, -1)))
census_1857$dist2 <-census_1857$dist1*(census_1857$krajna6_distance/1000)

census_1857$quad1<-ifelse((census_1857$lon > 18 & census_1857$lon< 20) & 
                            (census_1857$lat<47  & census_1857$lat>44), 1, 0)

census_1857$quad2<-ifelse((census_1857$lon > 16 & census_1857$lon< 18) & 
                            (census_1857$lat<47  & census_1857$lat>44), 1, 0)

census_1857$quad3<-ifelse((census_1857$lon > 14 & census_1857$lon< 16) & 
                            (census_1857$lat<47  & census_1857$lat>44), 1, 0)


census_1857$bfe1 <- ifelse(census_1857$krajna6_NEAR_FID == 1, 1,0)
census_1857$bfe2 <- ifelse(census_1857$krajna6_NEAR_FID == 2, 1,0)
census_1857$bfe3 <- ifelse(census_1857$krajna6_NEAR_FID == 3, 1,0)
census_1857$bfe4 <- ifelse(census_1857$krajna6_NEAR_FID == 4, 1,0)
census_1857$bfe5 <- ifelse(census_1857$krajna6_NEAR_FID == 5, 1,0)
census_1857$bfe6 <- ifelse(census_1857$krajna6_NEAR_FID == 6, 1,0)
census_1857$bfe7 <- ifelse(census_1857$krajna6_NEAR_FID == 7, 1,0)
census_1857$bfe8 <- ifelse(census_1857$krajna6_NEAR_FID == 8, 1,0)


census_1857_treat<-subset(census_1857, treat==1)
census_1857_treat2<-subset(census_1857, kreis_comitat_regiment_latin == "Militar-Grenz-Regiments-Bezirke" |
                             kreis_comitat_regiment_latin == "Militar-Grenz-Communitaten")


census_1857$population<-census_1857$catholic_latin+
  census_1857$catholic_greek+census_1857$catholic_armenian+
  census_1857$orthodox_greek+census_1857$orthodox_armenian+
  census_1857$evangelical_lutheran+
  census_1857$evangelical_reformed+census_1857$unitarian+census_1857$jewish

census_1857$pct_absent_1857<-census_1857$absent/census_1857$population*100
census_1857$POINT_X<-census_1857$lat
census_1857$POINT_Y<-census_1857$lon


pct_absent_1857<-lm(pct_absent_1857~treat + 
                      POINT_X + POINT_Y + 
                      zagreb_distance +
                      bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                      quad1 + quad2 + quad3, data = census_1857)
summary(pct_absent_1857)


###########
#2011 Data#
###########
data <- read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)
data$treat <- as.numeric( with (data,ifelse((data$treat_KRAJ_kraj6==1), 1, 0)))
data$dist1 <- as.numeric( with (data,ifelse(data$treat==1, 1, -1)))
data$krajna6_NEAR_FID
data$dist2 <-data$dist1*data$krajna6_distance


HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)
final_sp<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))


pct_absent_1944<-lm(pct_absent_1944~treat + 
                      POINT_X + POINT_Y + 
                      zagreb_distance +
                      bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                      quad1 + quad2 + quad3, data = data)
summary(pct_absent_1944)


migration_rate<-lm(migration_rate~treat + 
                      POINT_X + POINT_Y + 
                      zagreb_distance +
                      bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                      quad1 + quad2 + quad3, data = data)
summary(migration_rate)


######################
#Step1: RELEVANT VARS#
######################

relevant_vars<-c("pct_absent_1857",
                 "pct_absent_1944",
                 "migration_rate")

#####################
#Step2: MEANS AND SD#
#####################

mean_pct_absent_1857<-round(mean(census_1857$pct_absent_1857, na.rm=T), digits = 3)
mean_pct_absent_1944<-round(mean(data$pct_absent_1944, na.rm=T), digits = 3)
mean_migration_rate<-round(mean(data$migration_rate, na.rm=T), digits = 3)

sd_pct_absent_1857<-round(sd(census_1857$pct_absent_1857, na.rm=T), digits = 3)
sd_pct_absent_1944<-round(sd(data$pct_absent_1944, na.rm=T), digits = 3)
sd_migration_rate<-round(sd(data$migration_rate, na.rm=T), digits = 3)


means<-c(mean_pct_absent_1857,
         mean_pct_absent_1944,
         mean_migration_rate)

sds<-c(sd_pct_absent_1857,
       sd_pct_absent_1944,
       sd_pct_absent_1944)


####################
#Step3: LIST MODELS#
####################

list_models<-list(pct_absent_1857,
                  pct_absent_1944,
                  migration_rate)


######################
#Step5: Column labels#
######################


column_labels= c("\\shortstack{Pct. Out-Migration\\\\ 1857}",
                 "\\shortstack{Pct. Out-Migration\\\\ 1944}",
                 "\\shortstack{Net Migration Rate\\\\ 2011}")


croatia_grid2<-stargazer(list_models,
                        column.labels= column_labels,
                        keep = c("treat"),
                        covariate.labels=c("Military Territory"),
                        dep.var.caption = "Dependent variable",
                        dep.var.labels.include = F,
                        omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                        add.lines=list(c("Mean", means),
                                       c("SD", sds),
                                       c("Boundary FE", rep("Yes", 7)),
                                       c("Region FE", rep("Yes", 7))))
note_latex <- "\\multicolumn{4}{l} {\\parbox[t]{15cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and County clustered robust standard errors in parantheses from OLS regression. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The unit of analysis for the 1857 census can be city or district (stadt) or (land bezirke). 
The unit of analysis for the 1944 Catholic parish census is parish (zupa) and 
the unit of analysis for the 2011 Census is municipality (opcina).
See codebook in the appendix for data sources and how the variables were calculated.}}"
croatia_grid2[grepl("Note",croatia_grid2)] <- note_latex
cat(croatia_grid2, file="./Dropbox/Legacies_Project/Paper/tables/table_A4.tex", sep="\n")

